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ABSTRACT 

Models for the formation and growth of structure in a cold dark matter dominated 
universe predict that galaxy halos should contain significant substructure. Studies of 
the Milky Way, however, have yet to identify the expected few hundred sub-halos 
with masses greater than about 10^ M©. Here we propose a test for the presence of 
sub-halos in the halos of galaxies. We show that the structure of the tidal tails of 
ancient globular clusters is very sensitive to heating by repeated close encounters with 
the massive dark sub-halos. We discuss the detection of such an effect in the context 
of the next generation of astrometric missions, and conclude that it should be easily 
detectable with the GAIA dataset. The finding of a single extended cold stellar stream 
from a globular cluster would support alternative theories, such as self-interacting dark 
matter, that give rise to smoother halos. 



1 INTRODUCTION 



cUSTTIUr 



A genelic piedlcLiun uf Culd Daik MaLLei (CDM) 
ogy is the existence of a significant amount of substructure in 
gravitationally collapsed structures suc h as galaxy clusters 



galaxy halos and even d warf galaxies (Klypin et al. 1999 



Moore et al. 1999t 



2001). On the scale of large spiral galax- 
ies like the Milky Way, some 500 dense clumps are expected 
to orbit in the Halo. These structures have very dense cores, 
with a very steep radial profile of univers al shape (increas- 



mgasr 

[White 19971 iMoore et al. 



in the central regions, Navarro, Fronk & 



1999a), which renders them essen- 



tially impervious to Galactic tides. 

The current "dark matter crisis" stems from the differ- 
ence between these predictions and observations of low-mass 
galaxies. The first problem is that rotation curves of gas- 
rich low- mass disk galaxies are not as peaked as required 



by CDM. Van den Bosch et al. (2000) suggested that this 



disagreement could be due to an error in the analysis of HI 
observations, in which the beam-smearing effect of extant 
HI radial velocity curves had not been taken into account. 
However, this appears not to be the case, as recent (high 



spatial resolution) Ha rotation curves have shown ( 


Cote et 


al. 200C; de Blok et al. 2001 




Marchesini et al. 2001). 



reioninzation on t he early evolutio n of small structure s (Bui 
lock ct al. 2OOOI ; ^omerville 200l| ; [Tully et al. 200l[ ). This 
re-ionizing radiation sets in at a redshift ~ 6-10, and origi- 
nates from either the first stars or quasars. Indeed, the first 
direct e vidence for reioniza tion (at z ~ 6) has recently been 
found (Becker et al. 2001). The intergalactic and galactic 
mediums are ionized by this radiation, and any gas that is 
not in deep potential wells is lost from the protogalaxies. 
The limiting mass for maintaining a gas component appears 
to correspond to a circular velocity of ~ 30kms~^, well 
above the circular velocity of all the dwarf spheroidal galax- 
ies known. (However, it may be possible for low-mass sys- 
tems to maintain a gas-fraction if the gas managed to cool to 
dense molecular form before the epoch of reionization). This 
"squelching" of star-formation in low-mass galaxies neatly 
solves the CDM satellite over-production problem. 



A potential candidate for the missing dwarf systems 
that are supposed to inhabit the Galactic halo is the large 
population of high velocity clouds that pepper the sky (see 
Wakker 2001 for a recent review). With little or no stel- 
lar content, these are seen to possess velo cities incompatible 
with Galactic rotation models, although Blitz et al. (199£) 



A second fundamental problem also exists: the mass of 
the CDM substructures in a galaxy-sized halo is similar to 
that inferred for galaxy satellites, yet the number of CDM 
substructures in a galaxy-sized halo is much larger than the 



number of observed dwarf galaxies( Klypin et al. 1999 



Moore 



et al. 1999b). It is, of course, possible that these galaxies 
have a quite different non-cosmological origin (perhaps due 
to tidal interactions in the early universe). In either case 
the pro blem remains: where is the population of 500 CDM 



have shown that the kinematic characteristics of the veloci- 
ties of the clouds are compatible with them being members 
of an accreting population distributed within a megaparsec 
of the Milky Way. At these distances, the kinematic sig- 
natures of the clouds im ply a mass-to-light ratio of 10-50 
Burton (2000| ). Currently, their distances are ex- 



Braun 



tremely difficult to determine, with Complex A being the 
only system to which the distance, of between 4 and lOkpc, 
has been measured ( van Woerden et al. 199e|). With the 
promise of acc urate distances from Ha (Bland-Hawthorn & 



substructures orbiting the Milky Way and other large galax- 
ies? Recently, several studies have investigated the effect of 



Maloney 1999), a full appraisal of the nature of high velocity 



clouds and their relation to the formation of the Milky Way 
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can be made; but a preliminary analysis suggests that the 



clouds are scatte red in the Galactic halo (Bland-Hawthorn 



& Malcney 2001), and are not of cosmological origin. 

Nevertheless, CDM theory predicts that a large popu- 
lation of completely or almost completely dark "galaxies" 
should be present in the Milky Way. Here we propose an 
observational test of this scenario. In the early galaxy, many 
globular clusters are expected to have surrounded the Milky 
Way. Galactic tidal forces have been destroying these ob- 
jects ( [Gnedin fc Ostriker 1997 ), and of the initial popu- 
lation of a few hundred objects, only ~ 10 remain, and 



many o f these show signs of tidal disrupt ion ( Leon, Meylan 
& Com laes 20o"^ ; [Odenkirchen et al. 200l| ). 

In a static potential, as a low-mass stellar system, such 
as a globular cluster, loses progressively more and more mass 
to tidal forces, it develops long tidal tails that closely follow 
the orbit of the globular cluster. (In the limit of a zero mass 
system, the tidal stream exactly follows the progenitor or- 
bit). Here we explore how the presence of the population 
of dark matter clumps predicted by CDM alters the phase- 
space structure of a globular cluster tidal stream. For the 
special case when the global Galactic potential is nearly 
spherical, this corresponds to a broadening of the stream 
from a thin great-circle stream into a wid e band on the sky. 
In a companion paper (Ibata et al. 2002; Paper 2), we use 



these results to motivate the search for tidal streams in the 
2 Micron All Sky Survey (2MASS) dataset. 

Other ways to probe the dark matte r subs tructure hav e 
been proposed. Recently, [Metcalf (200l| ) and |Chiba (200"l| ) 
suggested that substructure in external galaxies may be 
detected through their gravitational lensing of background 
quasars, as a lumpy halo will produce image and brightness 
configurations different to a smooth matter halo. While such 
an approach may reveal any missing sub-halos, the method 
is subject to the vagaries of gravitational lens modeling. Dur - 
ing the final preparation of this work, Mayer et al. (2001) 



presented a study investigating the signatures that are im- 
printed by different dark matter models on tidal streams of 
dwarf spheroidal galaxies. They find that the structure and 



kinematics of the outer regions of dwarf galaxies can also 
allow one to ditterentiate between dark matter models. 



2 NUMERICAL SIMULATIONS 

Our aim is to quantify the difference in the structure of 
the tidal tail of an ancient globular cluster (or even an an- 
cient globular cluster remnant) if the Galactic halo has a 
smooth mass distribution compared to the case when the 
Halo also contains significant substructure. We have mod- 
elled the smooth components of the Gala xy as fixed po- 
tentials, u sing the Galactic mass model of 
ney (199^), which contains a disk, thick disk 



Dehnen & Bin- 



interstellar 

medium, bulge, spheroid and halo components. The p aram- 
eters of this model are detailed in Ibata et al. (2001). The 
halo was taken to have a mass normalization that gives a to- 
tal model circular velocity of Vc = 210 kms^^ at 50kpc. The 
mass fiattening of the potential was left as a free parameter 
gm- 

We adapted the fast parallel code PKDGRAV to include 
the forces due to these mass distributions by including a 
multipole expansion code kindly supplied by W. Dehnen. In 



all the following simulations we maintained accurate forces 
by using an opening angle of 0.75 and expanding the cell 
moments to hexadecapole order. Two body relaxation was 
suppressed by using a spline softening of 10 pc, such that 
the inter-particle forces are completely Newtonian at 20 pc. 
A variable time-step scheme is used based on the local ac- 
celeration, < r;-^/[$[/a, and de nsity. At < r\l y/G p, with 
the accuracy parameter r\ — 0.03 (Quinn et al. 1997). Typ- 



ically we used 100000 base steps to integrate 12 Gyr. With 
the variable time-steps, this was equivalent to taking 5 x 10^ 
time- steps. 

Th e initial glo bular cluster was modelled with a King 
model (King 1966), populated with 10* particles. The total 
mass of the globular cluster was taken to be 10® M0 , the 
model had concentration parameter c = 1.0 and central ve- 
locity dispersion a = 4.1kms~^, yielding a tidal radius of 
rt = 300 pc. This model is of approximately the same size 
and mass as Omega Centauri, but of lower concentration (by 
a factor of 1.7). Note that we are not aiming at presenting an 
accurate model of the disruption of a globular cluster — we 
do not have sufficient resolution for that — instead, here we 
seek to study the kinematic behaviour of the unbound tidal 
streams that emerge from disrupting low-mass systems. 

The initial position and systemic velocity of the King 
model were ra ndomly chosen from the halo component of a 
galaxy model (Boily, Kroupa & Penarrubia-Garrido 2001) 



that gives a similar rotation curve to the adopted Dehnen 
& Binney mass model. 

Integrating the globular cluster models for 10 Gyr in 
the smooth potential gives rise to very narrow tidal tails. 
An example of one of these models is shown in Figure 1, 
which displays the position, heliocentric radial velocity and 
distance distribution of the stream particles. At its narrow- 
est, the width of this stream is ~ 100 pc, similar in width 
to the tidal radius of the initial King model, while at its 
widest, the stream is ~ 3 kpc wide. The shape of the halo 
mass distribution for the simulation displayed in Figure 1 is 
spherical, qm ~ 1, which is motiv ated by recent analyses of 



the outer h alo of the Milky Way (Ibata et al. 2001 



Beers 2001) 



Chiba & 



For each one of these models, we also investigate the 
effect of the presence of a large number of moving substruc- 
tures in the halo. However, we first need a realistic way of 
modeling these substructures. According to Navarro, Frenk 
& White (1997) (hereafter NFW), and other authors, galaxy 
halos have a "universal" density profile p(r) oc l/(r(l -l-r)'^), 
which fits all scales currently probed by the numerical cos- 
mology simulations. The potential corresponding to this 
density profile is very simple: (;6(r) oc 1 — ln(l -I- r)/r, and 
the forces due to this potential would be straight-forward 
to add in to our N-body integrator. Unfortunately, however, 
this density profile is not well-behaved either at small or at 
large radii (the total mass diverges); to fix these problems, 
the force near the center of the sub-halo would need to be 
softened, and the density would have to be made to fall off 
to zero at some large radius (the tidal radius of the sub- 
structure would be a good choice). A simpler solution is to 
approximate the NFW profiles with a softened point-mass 
potential. Figure 2 compares the radial acceleration due to 
an NFW potential (dotted line), a Kepler potential (dashed 
line) , and a softened Kepler potential (solid line) with spline 
softening of 3.41 times the scale radius of the NFW pro- 
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file. Note that, at all radii, this softened point-mass potential 
underestimates the forces compared to the NFW potential. 

The substructures are added into the Galactic halo as 
softened point masses, with a distribution of sub-halo masses 
that follows a relation loRmfA^e ) = 3(1 — 'ivc/VGio bai) (our 



parameterization of Figure 2 of Moore et al. 2001). At the 
low-mass end of this distribution, there are 435 sub-halos of 
circular velocity Vc > 0.03VGiobai (M ~ lO'^Mo). The halo 
mass fraction in this lumpy component is relatively small, as 
it accounts for less than 10% of the mass of the smooth halo. 
The initial positions and velocities of the particles represent- 
ing th e sub-halos are also chosen from the halo co mponent 
of the Boily, Kroupa & Penarrubia-Garrido (2001) Galaxy 



model, but with a pericenter cutoff of 10 kpc, since dynami- 
cal friction will likely make the orbits of dark satellites with 
smaller pericenter radii decay quickly. 

The effect of these numerous dark matter clumps is 
marked: in Figure 3 we show the position, velocity and 
distance distribution of the tidal stream particles. In this 
simulation, only 1% of the mass of the halo is not in the 
smooth component, yet the tidal stream has been substan- 
tially heated and is substantially wider than in the smooth 
halo (the stream is so fluffy that it is diflicult to measure a 
width) . 

A better way to look at the simulations is in the space 
of the integrals of the motion. In the axisymmetric poten- 
tials considered here, the integrals are the total energy per 
unit mass E and the z component of angular momentum 
per unit mass Lz . It is also useful to inspect the total angu- 
lar momentum per unit mass, L, which is an approximate 
integral of the motion in the relatively spherical potentials 
we consider here. Figure 4 shows the relation between these 
three quantities. The top row of Figure 4 is derived from 
the simulation of the globular cluster model in the smooth 
spherical halo whose structure on the sky was previously 
displayed in Figure 1. After 12 Gyr the particles have re- 
mained localised in the space of the integrals of the motion. 
The bottom panel shows the same globular cluster model on 
the same orbit, but simulated in the halo model that con- 
tains the additional NFW sub-halos (corresponding to the 
simulation shown in Figure 3). In this lumpier potential, 
the globular cluster stream is comparatively much more dis- 
persed, especially in the parameter. The r.m.s. dispersion 
in Lz is a factor of 5 larger in the lumpy halo simulation. 

Inspecting the spatial distribution of particles can be 
useful for differentiating between a smooth and lumpy halo 
when the mass distribution is spherical (or close to spher- 
ical). However, in a flattened halo, differential precession 
mimics the effect of heating by CDM clumps, and it be- 
comes impossible to differentiate between the two causes of 
spatial dispersion of the stream. It is only when viewing the 
space of the integrals of the motion that the effect of pre- 
cession decouples from the stochastic heating by the CDM 
clumps. We have repeated the same experiments as above in 
a flattened halo, with = 0.7, a flattening value predicted 
by r iunicrical cosmology (pC atz 199 1; [Pubinski fc Carlbcrg 
1991|; [Vfarren et al. 1992 ; [Katz fc Guiiri 1991|; [Summers 1993r 



Dubinski 199^ The sub-halos are chosen from a halo distri- 



bution function which is also flattened to qm ~ 0.7. Figure 5 
shows the relation between the integrals of the motion E 
and Lz and the approximate integral L. Again, the dispers- 
ing effect of the sub-halos is very strong, with the r.m.s. 



dispersion in increasing by a factor of 9 with respect to 
the smooth halo case. 



2.1 Suite of simulations 

To obtain a statistical understanding of the change in 
Lz, we undertook a series of simulations of the globular 
cluster model with initial position and velocity randomly- 
drawn from the spherical halo distribution function model 
mentioned above. Ten simulations were undertaken in the 
Galaxy model with the smooth spherical halo. A further ten 
simulations probed the same models with the addition of the 
NFW sub-halos, for a total of 20 simulations with qm = 1.0. 
To probe the situation with a more flattened halo, we un- 
dertook a further 20 simulations with qm = 0.7. 

The result of these experiments is shown in Figure 6, 
where we display the r.m.s. dispersion in Lz for simulations 
in a smooth halo (filled circles) and in in the presence of 
halo substructures (star symbols). The top panel shows the 
case in a spherical halo {qm = 1.0), while the bottom panel 
is for the flattened halo {qm ~ 0.7) simulations. 

Both in a spherical and in a flattened halo, the pres- 
ence of CDM-like substructures increases substantially, on 
average, the RMS dispersion in Lz of the stream stars. On 
some orbits (such as that of our simulation ^9, the effect 
can be small, however. Our simulations suggest that glob- 
ular cluster streams with a very dispersed specific angular 
momentum, say ctLz > 200kpckms~^ Mq"'^ require that 
the halo possesses significant substructure. Smaller values, 
(jZ/z < 50kpckms~^ M0~^ say, would require a smooth 
halo. 



3 OBSERVATIONAL CONSIDERATIONS 

We next consider whether it will be possible to identify a 
10^ Mq stream orbiting the Halo at distances up to 100 kpc 
(m — M — 20) with photometry, radial velocity, parallax 
and pro per motion data from the future astrometric mission 
GAIA (Ferryman et al. 2001). The great advantage of the 
GAIA dataset is that gives access to the full 6-dimensional 
phase space information for many stars, and at least 4- 
dimensional phase-space information (for those stars that 
are too far to have measured parallaxes and radial veloci- 
ties). In addition, the photometric information (in 15 bands) 
can give age and/or metallicity estimates, providing a pow- 
erful discriminant between stellar populations. 

Developing a full model of the Galaxy as viewed by 
GAIA is not a trivial task and it is beyond the scope of this 
paper. Here we simply make an order of magnitude estimate 
of the contrast of plausible globular cluster streams over the 
background by estimating the relative phase space density 
of stream stars over the large-scale Galactic spheroid (or 
stellar halo) component. We estimate the phase-space size 
of the stream as follows. The ^ lOkms"^ radial velocity un- 
certainty of GAIA exceeds the intrinsic radial velocity dis- 
persion of the stream, and is about 10% of the halo radial 
velocity dispersion. The proper motion dispersion at 100 kpc 
corresponds to about ~ lOkms"^ or 10% of the halo disper- 
sion (two coordinates). The positional width of the streams 



are ^ 3 kpc, or about 1% of the width of the halo. With the 
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full 15 photometric bands, the distance uncertainty of pho- 
tometric parallax measurements is likely to be < 10%. Thus 
the phase-space volume occupied by the stream is likely to 
be about one millionth of that of the halo. Since the stellar 
halo has a mass of only ~ 10^ M0 , the contamination will 
amount to approximately one star in a thousand. With a 
limiting magnitude ofV ^ 21, the GAIA measurements will 
easily reach down to the horizontal branch at 100 kpc. Plac- 
ing all of the stream at a distance of d = 100 kpc, we expect 
~ 1600 stars for V > 21 per 10*' Mq of progenitor mass (as- 
suming the luminosity function of the Draco by Odenkirchen 
et al. 2001, corrected to a mass to light ratio of M/L = 3). 
For radial velocity measurements the limiting magnitude is 
V ~ 17, and the number of available stream stars from the 
10^ Mq population drops to ~ 180. The detection of these 
streams will be straightforward with GAIA, and the increase 
in the dispersion of due to any halo substructure will be 
easily resolved. 

With results expected towards the end of the second 
decade of this century, the GAIA mission is not for the 
impatient. So it is worthwhile to investigate whether any 
presently available datasets already contain sufficient infor- 
mation to uncover the predicted Galactic halo dark sub- 
structures. To this end we undertook an analysis of the 
2MASS dataset, which provides positions and infrared CCD 
photometry for sources brighter than K = 14.3. This gives 
access to a relatively restricted sample of Galactic halo M- 
giants, with no information on their kinematics, and only a 
very poor constraint on their distances. The search for the 
halo streams in the 2MASS dataset is presented in Paper 2; 
the analysis reveals a single very strong stream, that due to 
the ongoing disruption of the Sagittarius dwarf galaxy. 

Is it possible to use the Sagittarius stream to constrain 
the lumpiness of the halo? To investigate this issue we ran 
two simulations, in the same manner as above, one in a 
smooth halo, and the other with NFW clumps, using as 
an initial model f or the Sagittarius dwarf the model D2 of 
[bata et al. (2001 ) . Figure 7 shows the result of this experi- 
ment. This model is much more massive (M = 5 x 10* M©) 
than the globular cluster model simulated above, and is ini- 
tially much more extended (half-mass radius of 0.9 kpc), so 
it is naturally widely dispersed in E, Lz and L. Though the 
stream particles are heated and dispersed by the sub-halos, 
the effect is subtle, and will be very difficult to disentangle 
from the normal dynamical evolution of the dwarf galaxy. 
So low mass streams from globular clusters need to be found 
to act as probes of the dark matter. 



4 CONCLUSIONS 

We have shown that stellar streams of low mass systems such 
as disrupted globular clusters can be used to constrain the 
substructure of the dark halo. If current cosmological simula- 
tions correctly represent reality, ancient stellar streams from 
globular clusters should be substantially dispersed over the 
sky, even if on average, the Galactic halo potential is close 
to spherical. This effect is readily noticeable, as the streams 
will have a large dispersion in the z-component of specific 
angular momentum L^. For the particular halo studied in 
this contribution aL^ > 200kpckms~^ M©"^ implies the 
presence of substantial dark matter substructure. 



We have simplified the analysis in two major ways. 
Highly softened point-masses were used to model the halo 
substructures, which leads to an underestimate of the in- 
creased dispersion in Lz in the presence of the NFW clumps. 
Our other major assumption, that the Galactic potential is 
axisymmetric and does not evolve with time, will likely have 
the opposite effect, namely to underestimate the dispersion 
in the absence of the NFW clumps. Our study should be re- 
fined once the Galactic potential is better constrained (with 
data from the GAIA and SIM missions). 

This work has shown that an all-sky survey with high 
precision astrometry, such as will be obtained with the 
GAIA satellite, we will easily be able to detect globular clus- 
ter streams, and measure their kinematic properties, espe- 
cially the dispersion in L^. 

The heating due to the numerous encounters with the 
dark substructure must also affect the computed destruction 
rate of globular clusters and low-mass dwarf galaxies alike. It 
is possible that the lower limit to the mass of dSph galaxies 
~ 10^ Mq may be set by these disruptive encounters. 

Discovering and determining the thickness of ancient 
stellar streams will allow us to probe the invisible dark halo 
structures predicted by Cold Dark Matter cosmology. Find- 
ing a single narrow stream of width comparable to the tidal 
radius of a globular cluster would stand in strong contra- 
diction with this theory. Data from the forthcoming GAIA 
mission planned by the European Space Agency will provide 
the opportunity to undertake this study. 
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Figure 1. The sky distribution (in Galactic coordinates) of one of 
the globular cluster models after 12 Gyr. The Galactic potential 
has been modelled as a sum of a disk, thick disk, interstellar 
medium, bulge and spheroid components, plus a spherical dark 
halo. In this simulation all of these Galactic components have 
potentials that are smooth, axisymmetric, and static. The upper 
diagram shows the particle velocities and the lower diagram shows 
the particle distances. 



Figure 2. The acceleration as a function of radius due to an 
NFW halo (dotted line) a point-mass (dashed line) and a spline- 
softened point mass with Tg^ft = 3.41rs. By modeling the NFW 
sub-halos as softened point masses, we underestimate the forces 
(and hence the heating) they impart on other particles in the 
simulation. 
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Figure 3. As Figure 1, but this time tlie globular cluster has been 
simulated in a Galactic potential replacing a small mass fraction 
of the halo with 435 moving NFW sub-halos. 
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Figure 4. The relation between the integrals of motion for the 

simulation displayed in Figure 1 (top row) and the simulation 
displayed in Figure 3 (bottom row). The addition of a small frac- 
tion by mass in lumpy sub-halos has completely smeared out the 
stream. 
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Figure 5. As Figure 4, but in a flattened halo with Qm = 0.7. 
Evidently, even in a flattened halo, it is relatively straightforward 
to distinguish the effect of the sub-clumps. 
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Figure 6. The dispersion in the z-component of angular mo- 
mentum per unit mass (cLz) in several simulations. Filled circles 

denote the smooth halo simulations, while star-symbols show the 
results with NFW sub-halos. The top panel is for those simula- 
tions with a spherical mass distribution, while the bottom panel 
shows the results in a flattened halo qm = 0.7. Clearly, if the 
Galactic halo has substantial substructure as predicted by cosmo- 
logical simulations, globular cluster streams should have a wide 
dispersion in Lz- 
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Figure 7. As Figure 4, but for the case of a stream which has 
already been detected: that of the Sagittarius dwarf galaxy. Here, 
the difference between the smooth and lumpy halos is much less 
marked due to the fact that the initial model is much more mas- 
sive, and has naturally a wide spread in E, Lz and L. 



